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Abstract: The purpose of this paper is to study the problem of estimating a com- 
pactly supported function from noisy observations of its moments in an ill-posed 
inverse problem framework. We provide a statistical approach to the famous 
Hausdorff classical moment problem. We prove a lower bound on the rate of con- 
vergence of the mean integrated squared error and provide an estimator which 
attains minimax rate over the corresponding smoothness classes. 
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1. Introduction 

The classical moment problem can be stated as follows: it consists in getting 
some information about a distribution /i from the knowledge of its moments 
/ x^dfi{x). This problem has been largely investigated in many mathematical 
topics, among others, in operator theory, mathematical physics, inverse spectral 
theory, probability theory, inverse problems, numerical analysis but also in a 
wide variety of settings in physics and engineering such as quantum spin models, 
speech processing, radio astronomy (see for instance Kay and Marple (1981), 
Lang and McClellan (1983) . . . ). We may cite the classical and pioneer books in 
the field (see Akhiezer (1965), Shohat and Tamarkin (1943)) which put emphasis 
on the existence aspect of the solution and its uniqueness. According to the 
support of the distribution of interest, one may refer to one of the three types of 
classical moment problems: the Hamburger moment problem whose support of 
/i is the whole real line, the Stieljes problem on [0, -|-oo) and finally the Hausdorff 
problem on a bounded interval. In this paper, we shall focus on the last issue 
under the inverse problem angle. 
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The Hausdorff moment problem which dates back to 1921 (see Hausdorff 
(1921)) occupies a central place in the field of inverse problems and has been an 
object of great interest in literature since then. For instance, the particular case 
when only a finite number of moments are known which is called the Truncated 
Hausdorff moment problem, has recently aroused much attention (see Talenti 
(1987), Fasino and Inglese (1996), Inglese (1995), Tagliani (2002)). Another 
interest and aspect of the Hausdorff moment problem lies in its very close link 
to the inversion of the Laplace transform when this latter is given at equidistant 
points ( see for instance Brianzi and Frontini (1991), Dung, Huy, Quan and Trong 
(2006), TagUani (2001), Taghani and Velazquez (2003)). In fact, by a simple 
change of variable, the problem of Laplace transform inversion is equivalent to 
the Hausdorff moment problem. More recently, Ang, Gorenflo, Le and Trong 
(2002) have presented the Hausdorff moment problem under the angle of ill posed 
problems, in a sense that solutions do not depend continuously on the given data. 
Nonetheless, until now, as far as we know, the statistical approach which consists 
in assuming that the noise is stochastic has been very little put forward and rarely 
raised. 

We consider in this paper a statistical point of view of the Hausdorff moment 
problem. We aim at reconstructing an unknown function from noisy measure- 
ments of its moments on a symmetric bounded interval in a statistical inverse 
problem framework. In reality, it is barely impossible to measure moments with 
having any corruption. Without loss of generality we may and will suppose 
that [—a, a] = [—1,1]. In practise, observations of moments appear especially 
in quantum physics (see for instance Ash and Mc Donald (2003) and Mead and 
Papanicolaou (1983)), image analysis (see Teague (1980)), engineering mechanics 
(see Athanassoulis and Gavriliadis (2002)) and all the references therein). We 
shall give now some concrete examples of the application of the problem of mo- 
ments, in which either one can have access directly to moments to reconstruct 
the unknown function at stake or the moment problem appears in the line of 
reasoning. 

Examples. 1. In quantum physics, a big issue aims at reconstructing a positive 
density of states in harmonic solids from its moments on a finite interval. It is 
explained how to measure these moments (p 2410 and 2411 in Mead and Papan- 
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icolaou (1983)) and (p 3078 and 3079 in Gaspard and Cyrot-Lackmann (1973)). 

2. In the field of quantum gravity, the transition probabilities for a Markov 
chain related to the causal set approach to modeling discrete theories of quan- 
tum gravity satisfy a moment problem (see Ash and McDonald (2003)). One has 
to measure those probabilities. 

3. In the context of non classical moment problem, we may cite the following 
example (see Ang, Nhan and Thanh (1999)). It deals with the determination 
of the shape of an object in the interior of the Earth by gravimetric methods. 
The density of that object differs from the density of the surrounding medium. 
Assuming a flat earth model, the problem consists in finding a curve x — > cr{x) 
in the half plane < a{x) < H,0 < x < 1, a{x) satisfying a non linear integral 
equation of the first kind of the form 



where f{x) is a given function. The nonlinear integral equation can be approxi- 
mated by the following linear integral equation in ip: 



with (p{x) = H—a{x), M is large enough and x > 0. By taking x = 1, 2, . . . , n, . . . , 
we get the following equivalent moment problem: 



where /x„ = 27r/(— M — n), n = 1, 2, . . . 

In the first examples cited above, the authors aim at estimating an unknown 
density from measurements of its moments. The following paper presents the 
case of the reconstruction of an unknown function which could be in particular 
a density of probability, when one has its moments corrupted with some white 
noise. This approach constitutes one among other ill-posed inverse problems 
point of view. 
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Estimation in statistics using moments has already been put forward (see 
Mnatsakanov and Ryumgaart (2005)), but the approach there was based on em- 
pirical moments and empirical processes, the results were expressed in terms of 
weak convergence whereas our paper is built upon an ill-posed inverse approach 
with minimax results. 

The estimation procedure we use is based on the expansion of the unknown func- 
tion through the basis of Legendre polynomials and an orthogonal series method. 
We establish an upper bound and a lower bound on the estimation accuracy of 
the procedure showing that it is optimal in a minimax sense. We show that the 
achieved rate is only of logarithmic order. This fact has already been under- 
lined by Goldcnshluger and Spokoiny (2004). In their paper, Goldenshluger and 
Spokoiny (2004) tackled the problem of reconstructing a planar convex set from 
noisy geometric moments observations. They pointed out that in view of recon- 
structing a planar region from noisy measurements of moments, the upper bound 
was only in the order of logarithmic rate. The lower bound has not been proved. 
In a second part, they consider reconstruction from Legendre moments to get 
faster rates of convergence. Legendre moments can be observed in the context 
of shapes reconstruction. In our present work, instead of considering a planar 
region, we deal with functions belonging to a Sobolev scale and we stay focused 
on the classic moments with respect to the monomials x'^. Moreover, recently 
in the context of long-memory processes obtained by aggregation of independant 
parameter AR(1) processes and in view of estimating the density of the under- 
lying random parameter, Leipus, Oppenheim, Philippe and Viano (2006) had 
to deal with a problem of moments. They obtained very slow logarithmic rate 
but without showing that this could be the best possible. In a certain way, our 
minimax results provide a piece of answer. 

One might question this chronic slow rate which seems inherent to moment 
problems. In fact, the underlying problem lies in the non orthogonal nature of 
the monomials x^. They actually hamper the convergence rate to be improved 
for bringing a small amount of information. This remark is highlighted in our 
proof of the upper bound. 

This paper is organized as follows: in section 2 we introduce the model and 
the estimator of the unknown function and we finally state the two theorems. 
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Section 3 contains the proofs. The last section is an appendix in which we prove 
some useful inequalities about binomial coefficients. 

2. Statement of the problem 

2.1 The model. First of all, let us recall an usual statistical framework of 
ill-posed inverse problems (see Mathe and Pereverzev (2001)): 
Let A : H — > H he a known linear operator on a Hilbert space H. The problem 
is to estimate an unknown function f £ H from indirect observations 

Y = Af + e^, (1.1) 

where e is the amplitude of the noise. It is supposed to be a small positive 
parameter which tends to 0, ^ is assumed to be a zero-mean Gaussian random 
process indexed by on a probability space {Q,A, P). 

In the Hausdorff moment problem, the operator which determines the way 
the observations are indirect is defined by: 

Ak{f) = l'^x''f{x)dx, ken. (1.2) 

Let us now state the Hausdorff moment problem model. From the equation (jl.ip 
and ()1.2p . we derive the following sequence of moments observations pertubated 
by a stochastic noise: 

yk = P-k+£ik A; = 0,1,.... (1.3) 

where e is the noise level, it is supposed to tend to 0, are assumed to be i.i.d 
standard Gaussian random variables and are the moments of the unknown 
function / given by: 

= J x'^f{x)dx /c = 0, 1,.... 

The assumption that the noise is modelled as i.i.d standard Gaussian random 
variables is a natural modeling and has already been used (see Rodriguez and 
Seatzu (1993), Brianzi and Frontini (1991), Goldenshluger and Spokoiny (2004)). 
It is a standard assumption in statistical inverse problems. Moreover, we may 
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consider (|1.3p as a first approach of the problem at hand. 

The objective is to estimate the unknown function / supported on the bounded 
interval [—1, 1] from noisy observations of its moments in the model (II. 3p which 
can be assimilated to a gaussian white noise model but using a non orthonor- 
mal basis and which constitutes an inverse problem setting. The solution of this 
problem is unique in the case k = oo (see Bertero, De Mol and Pike (1985)), so 
that the statistical problem of recovering / from is relevant. In the partic- 
ular case of estimating a density of probability, the probability measure having 
density / is unique (see Feller (1968, Chap. 7)). 

The use of the Legendre polynomials in the HausdorfF classical moment prob- 
lem in order to approximate the unknown measure is a completely standard 
and quite natural procedure (see Ang, Gorenflo, Le and Trong (2002), Bertero, 
De Mol and Pike (1985), Papouhs (1956), Teague (1980), Goldenshluger and 
Spokoiny (2004)) as those polynomials directly result from the Gram-Schmidt 
orthonormalization of the family {x'^}, A; = 0, 1, . . . Consequently, expanding the 
function / to be estimated in the basis of Legendre polynomials falls naturally 
and fits the problem's nature. 

Denote Pnj the coefficients of the normalized Legendre polynomial of degree n: 

n 

Pn{x) — ^ ^ Pn,j 
3=0 

By considering the Legendre polynomials, if we multiply both sides of the model 
(II. 3p by the coefficients Pkj we get the following model (for the proof see Lemma 
[3l in Appendix): 

yk = Gk + s(^kCk (1-4) 

where al = Ej=o/?lj' Vk = J2'j=o PkjVj, Ok = J2j=oPkjf^j = j\ f{x)Pk{x)dx 
and are i.i.d standard Gaussian random variables. So the model (11. Sp is equiv- 
alent to the model (|1.4p . The model (|1.4p will be used for the proof of the lower 
bound. 

Before going any further, we can make a remark at this stage concerning the 
model (II. 4p which is an heteroscedastic gaussian sequence space model. Depend- 
ing on the asymptotic behavior of the intensity noise one may characterize 
the nature of the problem's ill-posedness (see Cavalier, Golubev, Lepski and Tsy- 
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bakov (2004)). Here, in our case, cr^ > j4l' 



(see Lemma [5l from Appendix) and 



hence tends to infinity exponentially. We may say that we are dealing with a 
severely ill-posed problem with log-rates . 

We assume that / belongs to the Sobolev space defined by: 



where 0^ = J^if{x)Pk{x)dx is the Legendre Fourier coefficient and P^, denotes 
the normalized Legendre polynomial of degree k. Note that we consider more 
general functions than densities of probability which are included in this smooth- 
ness class. 

Sobolev spaces associated with various kinds of underlying orthonormal basis 
constitute quite standard smoothness assumption classes in classical ill-posed 
problems (see for instance Mair and Ruymgaart (1996), Mathe and Pereverzev 
(2002), Goldenshluger and Pereverzev (2000)). In the Hausdorff moment prob- 
lem, the underlying basis is the Legendre polynomials. 

Let us now give some highlights of the Sobolev space VFJ regarding Legendre 
polynomials. Rafal'son (1968) and Tomin (1973) have shown in the more general 
case of Jacobi polynomials (and thus in the particular case of Legendre polyno- 
mials we are considering here) that Sobolev space W2 consists of all functions 
/ which have their derivatives /',/",.••> f^^~^^ being absolutely continuous on 
each interval [a,b] C (—1,1) (see also Mathe and Pereverzev (2002)). 

2.2 The estimation procedure. Let us define now the estimator of /. This 
latter is induced by an orthogonal series method through the Legendre polyno- 
mials. 

Any function in L^[— 1, 1] has an expansion: 



fc=o 

The problem of estimating / reduces to estimation of the sequence {0k}k=i 
Legendre polynomials form a complete orthogonal function system in L^[— 1, 1]. 



= if G L-'i-l,!] -.Y^k'^-m' < oc} 



k 
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We have 




Ok= f{x)Pk{x)dx = Yfik,j / f{x)x^dx = "^PkjfJ^j 



i=o -^-^ j=o 



This leads us to consider the following estimator of 6: 



k 




and hence the estimator /jy of /: 



N k 



fNix)=Y,OkPk{x) 



Y^PkjVjPkix) 



k=0 



k=0 j=0 



where yj is given by (jl.Sp and N is an integer to be properly selected later. 
The mean integrated square error of the estimator /jy is: 



where Ej- denotes the expectation w.r.t the distribution of the data in the model 
(jl.3p and for a function g £ L^[— 1, 1], 



In this paper we shall consider the problem of estimating / using the mean 
integrated square risk in the model (|1.3p . 

We state now the two results of the paper. The first theorem establishes an 
upper bound. 

Theorem 1. For a > 0, define the integer N = [a log (l/e)J ■ Then we have 



where C is an absolute positive constant and [-J denotes the floor function. 

We recall that the floor function is defined by: \x\ = max{n G Z| n < j;}. 
The second theorem provides a lower bound. 



%ll/jv-/f, 
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sup ¥.f\\fN- ft <C[\og{l/e)\ 



2r 
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Theorem 2. We have 

irif sup E;||/-/||2>c[log(l/e)]-2^ 

where c is a positive constant which depends only on r and the infimum is taken 
over all estimators f. 

3. Proofs 

3.1 Proof of Theorem [H For the following proof, we consider the genuine 
model (|1.3|) . By the usual MISE decomposition which involves a variance term 
and a bias term, we get 

N 

Ef\\fN-f\\' = EfY,i^,-9kf+ 

k=0 k>N+l 

but 

N N k 

EfY^{9k-ekf = EfY^{J2PkAyj-H)f 

k=0 k=0 j=0 

N k 

k=0 j=0 

and since ~ iV(0,l), it follows that 

N k 

k=0 j=0 k>N+l 

= Vn + B% 



We first deal with the variance term Vn- To this end, we have to upper bound the 
sum of the squared coefficients of the normalized Legendre polynomial of degree 
k. Set = Yl'j=oPkj- explicit form of Pk{x) is given by (see Abramowitz 
and Stegun (1970)): 



Pk{x) = 
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where [•] denotes the integer part and (^•) denotes the bmomial coefficient, 
(^) = 77-^- This involves 

[k/2] 



j=0 



< 



2k + 1 


1 


2 


4fc 


2k + 1 


1 


2 


4fc 


2k + 1 


1 


2 


4fc 



k\ f2k -2j 



2k\ 1 ^ r /k 



j=0 
\ I ,^k\2 



< ^^71^ ( J ^ i^r, (3.2) 



the inequality (|3.ip is due to the fact that for < j < [k/2], we have i^^ f^^'') < 
{^^). As for (1321), we have Efio' < { Ei=o (J)}' and it is well known 

that {E,to (?)}' = (2'=)^^- 

By using now that {(^^)} < ^ (see Lemma HI from Appendix) we have 



2 2A: + 1 
which yields 



Vn < Ce^N^lH^^, 

where C > denotes an absolute positive constant. 
Now, it remains to upper bound the bias term B'j^. 

k>N+l 



k>N+l 

00 

k=l 

Since the function / belongs to the space W2, Ylj^ ^^'^l^fcP < 00, we get 

Bj^ = 0(A^-2^) 

Finally we have the upper bound for the MISE: 
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At last, it remains to choose the optimal N which will minimize the expression 
(j3.3p . This is obtained by equalizing the upper bounds of the bias and the 
variance term, namely: 



as 4^^ ^> A^^'+'^/^, consequently x log(^). Once one plugs N x log(p-) in 
()3.3p . the desired result of the Theorem 1. follows. 

3.2 Proof of Theorem 2. From now on, to prove the lower bound and for prac- 
tical reasons, we shall consider the model (|1.4p which constitutes a heteroscedas- 
tic gaussian sequence space model. We recall that the equivalence between the 
models ()1.3p and ()1.4p is proved in Lemma El in Appendix. 
A successful approach and standard tool to obtain lower bounds for minimax risk 
consists in specifying a subproblem namely constructing a subset of functions 
based on the observations (jl.4p . Then we lean on the application of the following 
particular version of Fano's lemma (see Birge and Massart (2001)) which will 
allow us to evaluate the difficulty of the specified subproblem and will give us a 
lower bound for the MISE associated to this subproblem. 

One crucial point in the Fano's lemma is the use of the Kullback-Leibler diver- 
gence Ar(Pi,Po) between two probability distributions Pi and Pq defined by: 



Here's the version of Fano's lemma, we are going to exploit: 

Lemma 1. Let r] be a strictly positive real number, C be a finite set of elements 
{/o, . . . , fu} on M with \C\ > 6 and {Pj}j^c o. set of probability measures indexed 
by C such that : 



— 2r 




(i) U-fjW >??>0, VO<^<i<M. 

(ii) Pj <C Po, Vj = 1, . . . , M, and 



A'(Pj,Po) <H< logM 



then for any estimator f and any nondecreasing function i 
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First of all, we have to construct an appropriate set of functions £. We are going 
to define £ as a set of functions of the following type 

. 2m- 1 

k=m 

S={5^,...,S2m-i)eA = {0,iy 

We verify that fg belongs to • ^^^^ have to calculate the Legendre- 

Fourier coefficients associated with the function fg: 



051 = J Jg{x)Pi{x)da 



hence 



. ;(2.+2)/2 . 5^ if/eK2m-l] 
else 



+00 2 2m— 1 

k=0 k=m 
2 2m- 1 

k=m 



(3.4) 



k=m 

since 5^ G {0, 1}. 

We set (5^°) = (0, . . . , 0) and /^(o) = /o- The Legendre-Fourier coefficients of 
/o are null: 

Ooi = V/ G N. (3.5) 

We are now going to exhibit the suitable subset of functions C of the Lemma 1. 
To this purpose, we only take into consideration a subset of M + 1 functions of 
S: 

C = {fg(0) 5 • • • ) /5(*f) } 
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where {6^'^\ . . . , is a subset of {0, 1}'". 

We precise that is the law of the vector of observations Y = (yi, . . . ,yoo) in 
the model ^ for f = fs,6eC. 

We are now going to apply Lemma 1. We first check the condition (i), 
accordingly, we have to assess the distance 1 1/^(0 — f^u) |p. By the orthogonality 
of the system {Pk}k and thanks to Parseval equality we get, for < i < j < M, 



2m- 1 



k=m 

2 2m- 1 

k=m 

2 2m- 1 

- rn?r+i y°k Ofc 

k=m 



where •) is the Hamming distance. We are going to resort to the Varshamov- 
Gilbert bound which is stated in the following lemma to find a lower bound of 
the quantity p((5(*\ (J^^^): 

Lemma 2. (Varshamov-Gilbert bound, 1962). Fix m > 8. Then there exists a 
subset . . . , <5(*^)} of A such that M > T^l^ and 

p{5^i)^5ik)) VO < J < < M. 

8 

Moreover we can always take S^^^ = (0, . . . , 0). 

For a proof of this lemma see for instance Tsybakov (2004), p 89. 
Hence 

\\fsi^^-fswf>icl)/i8m'l^v' 



We are now going to check the condition (ii) in Lemma 1. and evaluate the 
Kullback-Leibler divergence. It is well known (see for instance Birge and Massart 
(2001) p 62) that for the Kullback-Leibler divergence in the case of a gaussian 
sequence space model we have 

KiFs,Fo) = 4f;^^l_^. (3.6) 
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Hence, by virtue of (|3.4p . (j3.5p and (|3.6p . the Kullback-Leibler divergence between 
the two probabihty measures Pq and of the observations in the model ([L 
associated respectively with functions / = /o anf f = fs for all 5 G C satisfies 

2 2r?i-l ,2r+2i:2 
- £2 ^4r+3 2^ 2 

but thanks to Lemma O (see Appendix) we have 

1 1 



af - 4^-1 
which implies 

K(P5,Po < -^-^5 ^-n > (5/ < TT-r, rn < 

One chooses m = log(p-) so that 1/4™ = e^, hence: 

i^(P5,Po) <c222-+4^ 

and since m < 8 log M/ log 2 (see Lemma 2) 

222r+7 

i^(P5,Po) < ^-^logM 
log 2 

Eventually one can choose cq small enough to have c = "^^^ < 1. 

Since now all the conditions of the Lemma 1. are fulfilled, we are in position 
to apply its result with the loss function ^{x) = and ry = (co)/(2\/2m^'). 
Therefore, we derive that whatever the estimator /, 



supE,[ii/-/in>^ 



2 c 

3 logM 



> 



but from above we had m = log(^) which gives the desired lower bound. 
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4. Conclusion 

The two theorems of this paper show that in the problem of estimating a 
function on a compact interval from noisy moments observations, the best rate of 
convergence one can achieve, supposing Sobolev scale smoothness and considering 
the mean integrated squared error, is only of logarithmic order. In a future work, 
one could try to generalise this result for loss and may obtain faster rate of 
convergence if one assumes a more restricted smoothness class involving super 
smooth functions. Besides, one may consider an heteroscedastic gaussian noise 
instead of a white noise model. 



5. Appendix 

Lemma 3. The models lil.3\) and are equivalent. 

Proof. We recall the model (II. 3p : 

Vj = /"i + ^^j = j f{x)x^dx + eij. 



We are going now to multiply both sides of (II. Sh by the coefficient j3jk of the 
Legendre polynomial: 

dLSl) <^ PkjVj = Pkj J f{x)x^dx + ePkj^j 

3=0 j=0 ''-^ j=0 

^ '^f^kjyj= f{x)^(3kjxHx + £^(3kjij 

3=0 ■^-1 3=0 3=0 

^ X ^''^y^ " / f'yx)Pkdx + e X ^^o^i 



Let us set = ^'j=oPkjij- Since are i.i.d standard Gaussian random vari- 
ables, the random variable ^k follows a normal law with zero mean and variance 
equal to X]j=o/^fcj- Hence (|1.3p is equivalent to: 

Vk = Ok + scFkik 
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where cr^ = Ej=o/5fc,j' Vk = T,j=o l^kjVj, Ok = Y.j=o (^kjf^j = f\ f{x)Pk{x)dx 
and are i.i.d standard Gaussian random variables. 

Lemma 4. For all n > 1 we have: 

Proof. Let us prove (15. ip by recursion on n. The inequahty is clearly true for 
n = 1. 

Suppose (jS.ip true for a certain n >1. 

2{n + l)\ _ f2n\2{2n + l) ^ 4*^ 2(2n + 1) 



n+1 J \n J n + 1 n^/^ n + 1 
by recursion hypothesis. It remains to prove that 



4" 2(2n + l) ^ 4"+^ 

^ — TT7I- (5-2) 



n 



1/4 n+1 1)1/4' 



(lOD ^ ^(^" + 1) < 4 

" ni/4(n+l) - (n + 1)1/4 



2 



n+l/2n + 1^4 ^ 

-^—( —y<2^ 

n n +1 ' 

which is true because we have (n+^)^/'^ < and (n+^)'^ < ?7-(n+l) 

since | < n^/2 + n/4. This completes the proof. 

Lemma 5. For all n > 1, we have: 

al > 4"-i (5.3) 

where an is defined in p.4p . 

Proof. Firstly, let us recall the value of the noise intensity cr^: 
2 _ 2n + 1 1 [ /n\ /2n - 2j^ ^ ^ 

n f2n^ ^ 
- 4^1 ^ 
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And so, in order to prove (j5.3p it remains to prove that 

> — 1= n > 1. 

n J 2y/n 

We again use a recursion on n. 

The inequality (j5.3p is clear for n = 1. We suppose the property true for a certain 
n > 1 and we shall prove it at the rank (n + l). 

'2(n + l)\ /2n\2(2n + l) 



n+l 7 \^ J n + l 

4" 2(2n + l) 



> 



n + l 



4"+i 

> (5.4) 

the inequality (j5.4p is true because it is equivalent to 4n^ + 4n + 1 > 4n^ + 4n 
what we always have. 
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